Tackling the Fermionic Sign Problem in the Auxiliary-Field 

Monte Carlo Method 



G. Stoitcheva, 1 W. E. Ormand, 1 D. Neuhauser, 2 and D. J. Dean 3 

1 Lawrence Livermore National Laboratory, 
P.O. Box 808, L-4U, Livermore, CA 94551, USA 
^ \ 2 Department of Chemistry and Biochemistry, 

^ ■ University of California, Los Angeles, CA 90024, USA 

3 Physics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA 

Abstract 

We explore a novel and straightforward solution to the sign problem that has plagued the 

Auxiliary-field Monte Carlo (AFMC) method applied to many-body systems for more than a 

decade. We present a solution to the sign problem that has plagued the Auxiliary-field Monte 

Carlo (AFMC) method for more than a decade and report a breakthrough where excellent agree- 

ment between AFMC and exact CI calculations for fully realistic nuclear applications is achieved. 

£> ' This result offers the capability, unmatched by other methods, to achieve exact solutions for large- 

ly ; 

scale quantum many-body systems. 
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The quantum many-body problem is the foundation of much of modern physics and chem- 
istry, and one of the great challenges in theoretical physics is to develop a fully microscopic 
solution that includes the full range of quantum correlations. Traditionally, configuration- 
interaction (CI) methods, such as the nuclear shell model, were a method of choice. Valence 
particles occupy a set of single-particle orbitals and influence each other through an effective 
interaction. In this valence space is a set of Nd basis states, (pi, used to construct the eigen- 
states of the Hamiltonian H. CI reduces to a matrix-diagonalization problem by finding 
the eigenvalues of the matrix Hij = (<j>i\H\<j>j). While powerful, CI is a brute-force method 
facing substantial computational challenges. Diagonalization methods for sparse matrices 
typically scale as Njj 25 and current limits are Njj ~ 10 9 ~ 10 . But, the number of basis states 
increases dramatically with particle number. For two-component systems, such as nuclei, 
the number of basis states increases as 




N D » s | " s |, (1) 

where is the number of proton (neutron) single-particle states in the configuration 

space, and N^ n > is the number of proton (neutron) valence particles. Consequently, for 
mid-mass nuclei, where N^ n ^ ~ 40 and ~ 20 would be typical, the matrix dimension 

would be of the order 10 20 , which to solve would require a computer 10 12 times more powerful 
than any available today. 

Monte Carlo methods offer an attractive alternative to CI as their computational effort 
scales more gently with particle number. Indeed, Monte Carlo methods have been applied 
to a wide variety of many-fermion problems in physics and chemistry; with applications in 
condensed matter, nuclear structure, and lattice quantum chromodynamics (see Ref. [lj]). 
Unfortunately, Monte Carlo methods applied to fermionic systems generally suffer from 
the well-known sign problem (where the sampling weight function is not positive definite), 
which substantially limits their efficacy Here, we will address the sign problem with the 
Auxiliary-Field Monte Carlo (AFMC) method |2| based on the Hubbard- Stratonovich (HS 
transformation 3(. AFMC has had wide applications; including the Hubbard model 4] 



and nuclei 



The full power of the AFMC method has not been realized because 



of the sign-problem, and past nuclear physics studies have been limited to using semi- 
realistic interactions with good sign or an extrapolation method based on breaking up the 
Hamiltonian into "good-sign" and "bad-sign" parts [7|. Here, we report a breakthrough 



solution to the sign problem based on shifting the contour as initially proposed in Ref. 8|. We 
show that the shifted-contour method is practical and extraordinarily effective in mitigating 
the sign problem for fully realistic Hamiltonians in nuclear systems. These results offer a 
clear pathway to perform large-scale many-body calculations that include the full range of 
quantum correlations. 

The AFMC method for generic rotationally Hamiltonians is outlined in detail in Ref. jf| , 
and here we present the central features germane to our solution to the sign problem. AFMC 
is based on the imaginary-time evolution operator e~" to either filter from an arbitrary 
trial wave function, (f> , the ground-state (GS) value for the operator Q via 

i o |e-^/ 2 fie-^/ 2 |0 o ) 



<fi) 



GS 



lim 



or to compute the thermal expectation value 
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where Tr^jv) denotes the Z-proton and iV-neutron projected trace. Eqs. ([5]) and ([3]) are 
distinct and complementary approaches. Along with GS properties, the thermal formalism 
permits us to calculate structure information, such as electro-weak transition strengths, at 
finite temperature and is the optimal procedure for computing the density of states [9j; a 



10|. On the 



critical ingredient needed to describe astrophysical nucleosynthesis processes 
other hand, the "zero-temperature" formalism, Eq. ([2]), is an efficient way to compute GS 
observables. 

Since any two-body Hamiltonian may be written in quadratic form as 



(4) 



where here we choose <d a to be a generalized one-body density operator, V a the strength of 
the two-body interaction, and £„ the single-particle energies, we simplify exp(—j3H) making 
use of the (HS) transformation [3J] 
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where s = ±1 if A > or ±i if A < and o is the associated auxiliary field. Setting 
A = — pV a , we have 
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where G(a) = exp(— 1^*1 a a) * s the Gaussian factor, the volume element is T>[B] = 
Tlad-Va^P I V a | /27r, and /i(cr) = ^ Q ( £ a + ■Sa^ v cr a )0 Q ,. Since in general the operators O a 
do not commute, we split e~^ H into N t time-slices, i.e., e~^ n = ^ e ~ A i 3H ^ N t calculations 
presented here are with A/3 = 1/32 MeV -1 ) and apply the HS transformation at each time 
slice. Eqs. ©-(IB]) can then be written as 

6 _ fV[3]W(ff)(Cl) g 

{) ~ !D[a]W{o) ■ {7) 
Defining the one-body imaginary-time propagator as 

U a = e -^h{3(N t )) . . . e -Af3h(a(l)) ^ 

we express the weight function W(a) and as 

Ti( Z n) &U a 

W(a) = G(a)Tr {Z}N) [U a ], (Q) 9 = ' ' n J . (9) 

The auxiliary fields a a are not just parameters introduced for numerical convenience, 
but have physical significance. Their presence in h(a) essentially constructs a constrained 
mean field. Indeed, the maximum of the weight function, W(a), corresponds to the Hartree 
mean- field solution (|, satisfying the self-consistent condition 

at IF = -s a sgn(V a )(Q a ) S M F . (10) 

The principal advantage of the AFMC is that overall the computational effort scales 
much more gently with particle number. For example, the number of proton (neutron) 
auxiliary fields is at most (iVf( n )) 2 x Nt- Thus, while for the case where iVf = iV" = 40 
and N$ = N™ = 20 conventional CI methods are confronted with matrices with dimension 
~ 10 20 , the number of AFMC fields with N t = 100 is 320,000. 

Given the large number of auxiliary fields, Eq. (CO) is evaluated using Monte Carlo meth- 
ods. Thus, 

(^>MC = ^£<^> (11) 
i 

where N is the number of samples (typically 4000), and <?j is distributed according to W{a). 
The uncertainty in the integral is then governed by the variance. Central to the Monte 
Carlo evaluation is that W(a) be positive definite. For rotationally invariant applications, 
the general conditions for which W(a) > was examined in Ref. 6|, and was found to be 



true only for a small class of semi-realistic interactions, such as pairing-plus-quadrupole, for 
even-particle systems. Without a positive-definite weight function, we can try to proceed 
by sampling with |W(ct)|. Eq. (TTTi) is modified by the presence of the "sign", $(<?j) = 
W(ai)/ |W(<7j)|, multiplying (Cl)^, and is normalized by the average sign ($). Fig. 1 
demonstrates how AFMC fails for general Hamiltonians. The figure shows the thermal 
energy as a function of (3 for 28 Mg within the s<i-shell using the realistic Hamiltonian of 
Wildenthal 1 11. The solid line shows the exact CI result, where all 28,503 eigenvalues 
were obtained with OXBASH [12]. The circles show the AFMC calculation while sampling 



|IF(cr)| with the Metropolis algorithm [13J]. In general, sampling with |W(<t)I breaks down 
for (3 > 0.4 MeV 
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FIG. 1: Thermal energy for the nucleus 28 Mg computed with the sd-shell Hamiltonian of Wilden- 
thal. The solid line shows the exact CI result obtained from all 28,503 shell-model eigenvalues. 
The (blue) circles show the AFMC result using Metropolis sampling on |W(<r)|. The (red) triangles 
show the results obtained using the shifted-contour method. 

To address the sign problem, we rewrite the two-body Hamiltonian as 

E V M = E V«(Q a - °*Y + V a (2Q a a a - al), (12) 

a a 

and apply the HS transformation to the quadratic (O — a) 2 terms, giving for e~ Al3H 



A/3h(a) 
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where now h(a) = J2 a [ £ a + V a (s a a a + a a )]Q a . With the shift a a , the maximum of the 
weight function is now 

o a = -s a sga(V a )((Q a )g - a a ). (14) 

Thus, if we choose a a = the maximum of the weight function occurs at a a — 0. The 

presence of a a in the exponential factors in Eq. (fT3|) is important. For V a < 0, W(a) is 
shifted to the origin. While for V a > 0, the overall maximum is shifted into the complex 
plane with the maximum along the real axis at a a = 0. Further, a static phase is introduced 
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FIG. 2: Thermal energy and the state density p(E), for the nucleus 27 Na computed with the 
sd-shell Hamiltonian of Wildenthal. The solid line shows the exact CI result obtained from the 
eigenvalues. The circles show the AFMC result using the shifted-contour method. 



We note that since the maximum of W(a) is centered about a a = we can Monte Carlo 
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sample the cx-fields with the overall Gaussian factor G(a). The advantage of sampling with 
the Gaussian factor is that it offers an efficient method to sample uncorrelated values In 
Fig. [TJ the results of AFMC calculation of the thermal energy for 28 Mg using the shifted- 
contour method with Gaussian sampling is shown (triangles) and compared to the exact CI 
result as well as with Metropolis sampling on |VK(<t)|. Shifting the contour yields agreement 
with the exact thermal calculation, which clearly represents a significant improvement over 
previous capability. With the zero-temperature formalism, at (3 = 3.0 MeV -1 we compute a 
GS energy of -120.370(25) MeV, which is in good agreement with the CI result of -120.532 
MeV. 

In Fig. [2J we compare the CI and AFMC results for the thermal energy and the state 
density, p(E), for 27 Na in the sci-shell using the Wildenthal interaction. With the zero- 
temperature formalism, we obtain the GS energy of -99.106(55) MeV, which is also in good 
agreement with the exact result of -99.230 MeV. The state density (the total density of 
states including the (2J + 1) degeneracy for each state of angular momentum J) can be 
computed with the saddle-point approximation for the inverse Laplace transform of the 
partition function, i.e., 

p(E) = e ln z^+^W /^/-2ndE((3) /dp, (15) 

where In Z(/3) = — Jjf d/3'E(j3') + lnZ(0), and Z(0) is the total number of states given by 
Eq. dU). These 27 Na results are significant because not only is excellent agreement with the 
CI calculation achieved, but beforehand odd-systems suffered from the sign problem even 

n 

for the semi-realistic, good-sign interactions identified in Ref. [6J. 

In Fig. [3l w e show results for the more challenging case of 56 Fe, where the GXPF1A 
interaction [14[] was used in an active model space comprised of the 0/ — lp orbits. Here, 
N p(u) = 20; np = 6, and N" = 14, and the number of CI basis states with J z = is « 501M. 
In this case, we obtained the GS energy with the shell- model code REDSTICK [ijj], which is 
represented by the solid line in the figure. The shifted-contour AFMC calculation is clearly 
converging to the full-space CI result. With the zero-temperature formalism we calculate 
a GS energy of -195.687(107) MeV, which is in good agreement with the CI result of - 
195.901 MeV. The computational advantage of AFMC for large model spaces is evident as 
the zero-temperature calculation for 56 Fe took 12 CPU hours [16], as opposed to 1000 CPU 



hours for CI 



17J]. In the upper panel, we compare the calculated state density with values 



i 1 1 r 




a Ref. [19] 
Ref. [18] 



E (MeV) 



CO. 



■140 



■160 



■180 



i i i i i i i i i i i i i i i i i i i i i i i i 



• AFMC - Shifted Contour 
CI Ground state 




_9QQ I i i i i I i i i i I i i i i I i i i i I i i i i 

0.5 1 1.5 2 2.5 
P (MeV" 1 ) 

FIG. 3: Thermal energy and the state density p{E) for the nucleus 56 Fe computed with the 
GXPF1A /p-shell Hamiltonian. The solid line in the bottom panel shows the exact CI result for 
GS energy. The circles show the AFMC result using the shifted-contour method. In the upper 



panel, the calculated state is compared with values inferred from recent experiments [(squares) 
(triangles) Q]. 



inferred in recent experiments 
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19]. Overall agreement with the inferred experimental 
quantities is achieved. Here, our intent is to demonstrate a new capability, thus our AFMC 
calculation consists of just one major shell. Consequently, negative-parity and higher-lying 
states are outside this model space, and the calculated state density will under predict the 
observed state density at higher excitation energies. In principal, there are no underlying 
computational difficulties in extending our calculations to include more major shells; only 
the question of the appropriate effective interaction. 
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We present a solution to the sign problem for the AFMC method applied to many-body 
systems based on shifting the quadratic part of the two-body Hamiltonian. The optimal 
choice for the shift is the fields associated with the Hartree mean-field solution for each 
specific value of (3. This choice shifts the maximum of the integrand to the origin; permitting 
efficient sampling using the Gaussian factor. For bad sign components of the Hamiltonian, 
the shift introduces phases that mitigate the presence of negative signs in the weight function 
as the fields are sampled along the real axis. With A/3 = 1/32 MeV -1 , the thermal energy is 
typically reproduced at the level of 300 keV or better, while the GS energies are reproduced 
to within 150-200 keV. This is a substantial improvement over previous attempts [7|, where 
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ill, Q. 



deviations of the order 1 MeV from CI results were common 
generally the level of accuracy achieved by the effective interactions themselves 
An exciting possibility for the future is to combine the AFMC with traditional mean-field 
approaches based on Skyrme-like interactions [2l| to develop a universal picture for nuclei. 
In addition, given that the mean field is ubiquitous for many-body systems, it is likely that 
the shifted-contour method will have wide ranging applications to the quantum many-body 
problem across many subfields of theoretical physics and chemistry. 
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